Chapter 8 Beta diversity

load("data/data.Rdata")
beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),] 
rownames(genome_counts_filt) <- NULL
genome_gifts <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]

dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

8.1 Permanova

#Richness
betadisper(beta_q0n$S, sample_metadata$region) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.09234 0.092344 14.312    999  0.001 ***
Residuals 56 0.36132 0.006452                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Daneborg Ittoqqortoormii
Daneborg                             0.001
Ittoqqortoormii 0.00037875                
adonis2(beta_q0n$S ~ region,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_w3naam8l0g4kwhl0jbuc
term df SumOfSqs R2 statistic p.value
region 1 1.700798 0.2101443 14.89903 0.001
Residual 56 6.392678 0.7898557 NA NA
Total 57 8.093476 1.0000000 NA NA
#Neutral diversity
betadisper(beta_q1n$S, sample_metadata$region) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq     F N.Perm Pr(>F)    
Groups     1 0.39850 0.39850 26.96    999  0.001 ***
Residuals 56 0.82775 0.01478                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Daneborg Ittoqqortoormii
Daneborg                             0.001
Ittoqqortoormii 2.9964e-06                
adonis2(beta_q1n$S ~ region, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_5v9x7py4dufh29pjo842
term df SumOfSqs R2 statistic p.value
region 1 2.376111 0.2881769 22.67123 0.001
Residual 56 5.869211 0.7118231 NA NA
Total 57 8.245323 1.0000000 NA NA
#Phylogenetic diversity
betadisper(beta_q1p$S, sample_metadata$region) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.10517 0.105171 12.067    999  0.001 ***
Residuals 56 0.48808 0.008716                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Daneborg Ittoqqortoormii
Daneborg                             0.002
Ittoqqortoormii 0.00099752                
adonis2(beta_q1p$S ~ region, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_cu9ouyop6wizsg7srdqf
term df SumOfSqs R2 statistic p.value
region 1 0.1387067 0.101621 6.334495 0.001
Residual 56 1.2262339 0.898379 NA NA
Total 57 1.3649405 1.000000 NA NA
#Functional diversity
betadisper(beta_q1f$S, sample_metadata$region) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003452 0.0034517 0.8969    999  0.353
Residuals 56 0.215512 0.0038484                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                Daneborg Ittoqqortoormii
Daneborg                           0.355
Ittoqqortoormii  0.34768                
adonis2(beta_q1f$S ~ region, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
tinytable_ch63d6uhtyo5600panu4
term df SumOfSqs R2 statistic p.value
region 1 0.01958829 0.05063651 2.98689 0.177
Residual 56 0.36725289 0.94936349 NA NA
Total 57 0.38684118 1.00000000 NA NA

8.2 Plot

8.2.1 Richness diversity

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(region) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
      scale_color_manual(name="Region",
                         values=c("#6A9AC3","#F3B942")) +
      scale_fill_manual(name="Region",
                        values=c("#6A9AC350","#F3B94250")) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )

8.2.2 Neutral diversity

beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(region) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
      scale_color_manual(name="Region",
                         values=c("#6A9AC3","#F3B942")) +
      scale_fill_manual(name="Region",
                        values=c("#6A9AC350","#F3B94250")) +
    geom_point(size = 4) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )

8.2.3 Phylogenetic diversity

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(region) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
      scale_color_manual(name="Region",
                         values=c("#6A9AC3","#F3B942")) +
      scale_fill_manual(name="Region",
                        values=c("#6A9AC350","#F3B94250")) +
    geom_point(size = 4) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )

8.2.4 Functional diversity

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(region) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
      scale_color_manual(name="Region",
                         values=c("#6A9AC3","#F3B942")) +
      scale_fill_manual(name="Region",
                        values=c("#6A9AC350","#F3B94250")) +
    geom_point(size = 4) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )